function [model,kbasis] = cgparse6S_2(varargin)

% CGPARSE0 Solve for PARSE parameters using a variant of CG
%
% CGPARSE0(FID,CGP,MODEL)
%
% Pass data parameters in with the structure FID, CG algorithm control
% parameters with the structure CGP, and initial model parameter estimates
% with MODEL.  Output is returned in a model structure.  If coded correctly,
% the model should be determined from the initial estimate structure (future work).
%
% CGPARSE0(FID,CGP,MODEL,KBASIS)
%
% Also pass in precomputed KBASIS structure.
%

%------------- start of initialization section -------------------------
% input arguments
fid = varargin{1};
cgp = varargin{2};
model = varargin{3};
% pseq = varargin{4};
Ncpus = cgp.Ncpus;

if cgp.time
    tic;
end
if (nargin==4)
    kbasis = make_kbasis(fid,model,cgp.mat);
else
    kbasis = varargin{5};
end

pv = model.pvec;
CC = zeros(1,1000);  % set up variables for performance tracking
amp = zeros(model.reso,model.reso);expo1=amp;

% initialize the parameter array values as complex array
pv((1:model.N)+2*model.N) = single(real(pv((1:model.N)+model.N))/model.wfr + 1i*imag(pv((1:model.N)+model.N))/model.wfi);
%pvm=CC;pvr=CC;pvi=CC;
NF0 = 0;
% NM=model.N;
niter = 1;
CCold = 1e+15;
CCnew = 1e+14;
nitermax=5000;
modelcurr = model;

%kacq=max(sqrt(fid.kx.^2+fid.ky.^2));
%kr = sqrt(fid.kx.^2+fid.ky.^2);
%kf = max(kr);
kmsk1=single((fid.estwt));%windowing for signal model
%kmsk1 = (sqrt(1-(kr/kf).^2).*(abs(fid.kx)<kfs).*(abs(fid.ky)<kfs))';
%kmsk1=((abs(fid.kx)<kfs).*(abs(fid.ky)<kfs).*(0.5*(1+cos(pi*fid.kx/kfs))).*(0.5*(1+cos(pi*fid.ky/kfs))))';%...
%    .*(0.54+0.46*cos(pi*fid.kx/(48/22))).*(0.54+0.46*cos(pi*fid.ky/(48/22))))';
%[KX,KY] = meshgrid(((1:model.reso)-model.reso/2-1)/pseq.fovcm);
%kmask = sqrt(KX.^2+KY.^2) < kacq;

nloop = 0;
if cgp.time, toc, fprintf('end of initialization\n'), end
%------------- end of initialization section ---------------------------

if mod(niter,cgp.val) == 0,
    fprintf('%d %d %d %11.4e %11.4e\n',niter,cgp.NLIST(1),NF0,abs(CCold-CCnew),CCnew)
end
if mod(niter,cgp.time) == 0, toc, end

% begin iteration
tic;
for nlen = 1:cgp.NLENGTHS  % nested loops for progressive-length CG implementation
    fid.N = cgp.NLIST(nlen);  % current length
    fid.N = fid.N - mod(fid.N,Ncpus);
    modelcurr.pvec = pv;
    fid.est1 = kmsk1(1:fid.N).*fid1syn(modelcurr,kbasis,fid.N);%uses kbasis
    fid.est2 = kmsk1(1:fid.N).*fid2syn(modelcurr,kbasis,fid.N);%uses kbasis
    %fid.est = fidsyn(modelcurr,kbasis,fid.N);%uses kbasis

    while (abs(CCold-CCnew) > cgp.FLIST(nlen)*CCold || fid.N > NF0) && niter < nitermax
        nloop = nloop+1;  % keep track of iterations in each sub-loop
        gradJ = s_grad_2(modelcurr,fid,kbasis,cgp);
        if mod(niter,cgp.time) == 0, toc, end
        % define direction of 1-D search for SD or CG (CG = 1);
        if fid.N ~= NF0
            g = -gradJ;
%             gold = g;
            h = g;
        else
            gold = g;
            g = -gradJ;
            ga = max(0,real(((g-gold)*g')/(gold*gold')));
            h = g + ga*h;
        end
        %         if mod(niter,cgp.vis) == 0
        %             figure(13);title('gradients computed');
        %             plot([real(h); imag(h)].');drawnow;
        %         end

        modelcurr.pvec = pv;
        modelcurr.pvectmp = h;

        delpf = Newton2wtV_2(modelcurr,fid,kbasis,cgp);

        pv = pv + delpf*modelcurr.pvectmp;

        %         % bound frequency, and decay <=0
        if isfield(cgp,'bounddecay')
            ind2 = 1+2*model.N:3*model.N;
            binz = false(1,model.N);
            pvBndDecay = abs(model.wfr*real(pv(ind2))) > cgp.bounddecay;
            pv([pvBndDecay pvBndDecay binz]) = cgp.magdown*pv([pvBndDecay pvBndDecay binz]);
            pv([binz binz pvBndDecay]) = 0;
            pv(2*model.N+1:3*model.N) = min(real(pv(2*model.N+1:3*model.N)),0) + 1i*imag(pv(2*model.N+1:3*model.N));
        end
        modelcurr.pvec = pv;
        fid.est1 = kmsk1(1:fid.N).*fid1syn(modelcurr,kbasis,fid.N);
        fid.est2 = kmsk1(1:fid.N).*fid2syn(modelcurr,kbasis,fid.N);
        CC(niter) = jeval_2(fid,modelcurr,cgp);

        if mod(niter,cgp.plot) == 0
            figure(11);
            semilogy((CC(1:niter)));
            title(strcat(num2str(fid.N),'  --  ',num2str(CC(niter)),'    '));
            drawnow
        end
        for nl = 1:model.N
            amp(model.ix(nl),model.iy(nl)) = pv(nl);
            expo1(model.ix(nl),model.iy(nl)) = pv(nl+2*model.N);
        end

        R2e=real(-model.wfr*expo1)/fid.delt;
        frmap=imag(model.wfi*expo1)/(2*pi*fid.delt);
        if mod(niter,cgp.vis) == 0
            figure(14);
            subplot(2,2,1);
            imagesc(frmap);axis image
            title(strcat('deltaCC/CC = ',num2str(abs(CCold-CCnew)/CCnew)));
            subplot(2,2,2);
            imagesc(abs(amp));colormap hot;axis image
            subplot(2,2,3);
            imagesc(R2e);colormap hot;axis image
            title(strcat('est. R2*, ',num2str(fid.N)));
            subplot(2,2,4); plot(R2e(:,30));grid on        
            title(strcat('R2e profile ',num2str(niter),' iterations'));
            drawnow;
        end

        NF0 = fid.N;
        CCold = CCnew;
        CCnew = CC(niter);
        niter = niter + 1;
        if mod(niter,cgp.val) == 0
            fprintf('%d %d %d %11.4e %11.4e\n',niter,fid.N,NF0,(CCold-CCnew),CCnew)
        end

    end
    toc
    whos
end
toc
model.pvec=pv;

% varargout{1} = model;

if nargout==2
%     varargout{2} = kbasis;
end
if cgp.time, toc, fprintf('end of CG algorithm\n'), end
